rm(list = ls())

library(gsynth)
library(haven)
library(did)
library(panelView)
require(Rcpp) 
require(ggplot2)  
require(GGally) 
require(foreach)  
require(doParallel) 
require(abind) 
library(PanelMatch)
library(dotwhisker)
library(broom)
library(dplyr)

setwd("~/Dropbox/Immigration/Analysis/Gravity Models/")

data <- read_dta("gscm.dta")
options(scipen=999)

### SYNTHETIC CONTROLS: ALL ARRIVALS

policy <- gsynth(combinedarrivals ~ policytop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
policy

access <- gsynth(combinedarrivals ~ accesstop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
access

services <- gsynth(combinedarrivals ~ servicestop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
services

livelihoods <- gsynth(combinedarrivals ~ livehoodtop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
livelihoods

movement <- gsynth(combinedarrivals ~ movetop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
movement

participation <- gsynth(combinedarrivals ~ participtop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
participation

plot(policy, type = "raw", theme.bw = T, main="Policy Liberalization",
     ylab="FDP Arrivals", xlab="\n Year", xlim = c(2000, 2017), ylim=c(0, 800000))+
        scale_y_continuous(breaks=seq(0,800000,100000))+scale_x_continuous(breaks=seq(2000,2017,2))

plot(policy, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
     ggtitle("Policy Liberalization", subtitle="ATT: 35364, 95% CI: (18045, 52278)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(access, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Access Liberalization", subtitle="ATT: 29209, 95% CI: (9513, 48108)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(services, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Services Liberalization", subtitle="ATT: 29965, 95% CI: (9757, 52995)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(livelihoods, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Livelihoods Liberalization", subtitle="ATT: 28972, 95% CI: (7958, 49567)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(movement, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Movement Liberalization", subtitle="ATT: 31857, 95% CI: (14597, 50732)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(participation, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Participation Liberalization", subtitle="ATT: 29161, 95% CI: (13019, 46502)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))


### SYNTHETIC CONTROLS: TEK ARRIVALS

tekpolicy <- gsynth(tek_combinedarrivals ~ policytop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
tekpolicy

tekaccess <- gsynth(tek_combinedarrivals ~ accesstop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
tekaccess

tekservices <- gsynth(tek_combinedarrivals ~ servicestop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
tekservices

teklivelihoods <- gsynth(tek_combinedarrivals ~ livehoodtop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
teklivelihoods

tekmovement <- gsynth(tek_combinedarrivals ~ movetop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
tekmovement

tekparticipation <- gsynth(tek_combinedarrivals ~ participtop25_lag1+ ihs_gdppc_lag1+ polyarchy_lag1+ ihs_pop_lag1+ fariss_repression_lag1 +ihs_unemployilo_lag1 +cwbroad_lag1+ majoranncwbroad_reg1500_lag1, data = data, index = c("ccode2","year"), force = "two-way", CV = TRUE, min.T0=4, na.rm=T, r = c(0, 3), se = TRUE, inference = "parametric", nboots = 1000, parallel = FALSE, seed = 02139, conf.lvl = 0.95, EM = TRUE)
tekparticipation

plot(tekpolicy, type = "raw", theme.bw = T, main="Policy Liberalization",
     ylab="FDP Arrivals from TEK Origins", xlab="\n Year", xlim = c(2000, 2017), ylim=c(0, 800000))+
        scale_y_continuous(breaks=seq(0,800000,100000))+scale_x_continuous(breaks=seq(2000,2017,2))

plot(tekpolicy, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals from TEK Origins", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Policy Liberalization", subtitle="ATT: 36204, 95% CI: (15755, 55123)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(tekaccess, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals from TEK Origins", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Access Liberalization", subtitle="ATT: 34057, 95% CI: (14391, 55169)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(tekservices, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals from TEK Origins", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Services Liberalization", subtitle="ATT: 29808, 95% CI: (7663, 55202)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(teklivelihoods, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals from TEK Origins", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Livelihoods Liberalization", subtitle="ATT: 19921, 95% CI: (-3353, 41028)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(tekmovement, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals from TEK Origins", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Movement Liberalization", subtitle="ATT: 38062, 95% CI: (16123, 59018)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

plot(tekparticipation, type = "counterfactual", shade.post = T, theme.bw = T, main="",
     ylab="FDP Arrivals from TEK Origins", xlab="Time Relative to Treatment", xlim = c(-12, 12), ylim=c(0,250000))+
        ggtitle("Participation Liberalization", subtitle="ATT: 28807, 95% CI: (11339, 47045)")+theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))

